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Abstract 


A particle simulation study of the propagation and absorption of electro- 
magnetic waves and the associated plasma heating is presented. A plasma slab 
with nonuniform density profile, with positive gradient of density along the x 
-axis, is modeled with a two and one-half dimensional electrostatic particle 
simulation code. A finite length electrostatic single wave exciter is considered. 
The exciter frequency is chosen to be higher than but close to the lower hybrid 
frequency. Deceleration of the incident electromagnetic wave along the mag- 
netic field is necessary for its propagation. As the wave penetrates the plasma 
slab, its phase velocity reduces and at certain plasma density it transforms 
into a plasma wave. The conditions of this transformation, which have been 
derived theoretically using the hydrodynamic and the geometric o{)tics approx- 
imation are demonstrated computationally. The absorption of the waves by 
both ions and electrons is observed. Surface electron heating is also apparent 
in the simulations presented. Small sheaths are found near the exciter. .A.h 
predicted by the cold plasma theory, the resonance cones are observed in tin; 
phase-space plots. 
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N omenclat ure 


A list of symbols is given below. Some symbols which are specific to a section and whose 
meaning may be different in different sections are not included. Standard mathematical 
symbols also do not appear in this list. 
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Operation 


Vector derivative 

Latin Alphabet 


Bo 
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c 

Velocity of light 

D 

Displacement vector 

E 

Electric field strength 
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Zero*^ order velocity distribution for particle of type h 
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Reduced velocity distribution 
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J 

Current density 

k 

Wave propagation vector 

K 

Dielectric tensor 
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n 
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Wo 

x,y,z 
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Shape factor 
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Chapter 1 
Introduction 


The failure of fission to provide an environmentally clean energy source has 
motivated considerable interest towards nuclear fusion. The available informa- 
tion suggests that fusion energy can be a practically inexhaustible and cheap 
source of energy. After a long research that took place in various organizations 
worldwide, it was concluded that the nuclei of light elements can be fused ef- 
ficiently if the plasma of such elements is confined in a magnetic chamber and 
heated to its ignition temperature. Considerable amount of effort has been 
directed in the research of the magnetic confinement devices and the heating 
of plasma in such devices. Although the physics of plasmas was known well 
before the study of fusion reactions began, a different approach was required to 
attack the problems encountered in the development of fusion reactors. This 
lead to the evolution of a new branch of physics the Physics of Fusion Plasmas. 
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1.1 Role of Numerical Simulation in Under- 
standing the Physics of Fusion Plasmas 

At present plasma physics widely employs the numerical simulation techniques, 
particularly in the area of Heating. Plasma physics is based on sound math- 
ematical foundation. The equations describing the dynamics of the plasma 
were known before the fusion research began. However, fusion applications oi’ 
a plasmas require an analytic investigation of a mixed system, combining the 
equations of hydrodynamics, the kinetic equations, the Maxwell’s equations, 
and other auxiliary equations. The solutions of such a complex system can 
be approached analytically provided these equations are simplified to a large 
extent. On the other hand, the modern level of computcu's and achi('v('ments 
in numerical methods make it possible to analyze the given problem more n‘- 
alistically, considering the structural features of the exi)erimental (k'vice. th(‘ 
complex geometry of the magnetic field etc. Computers can be used to can'}' 
out the detailed quantitative analysis of a process in order to estimate the role 
of specific factors, reveal new properties and regularities, and considcu' comi)lex 
cooperative phenomena. 

The numerical methods are very important in comparing the theoretical and 
experimental results. A quantitative relation between them is difficult to estab- 
lish, as many quantities appearing in the theoretical model/results may not be 
observed or measured directly. The following lines are quoted by Dnestrovskii 

[Dnestrovskii and Kostomarov{\^i2)], 

« 

The theory and experiment use different languages and need an interpreter to 
communicate: this role is played by computers. 


Another role played by the computer simulation in fusion research is to make 
predictions. If a model has a reliable theoretical basis and is successfully tested 
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in practice, the accuracy of the predictions can be very high. Such predictions 
then become very important in planning the experiments. 

The complex nature of the problems that are encountered in the study of fusion 
plasmas has motivated large amount of interest in computer simulation, which 
has played an essential role in the development of the theory of fusion plasma. 
Computer simulation of plasmas comprises two general areas, based on either 
the kinetic description or the fluid description Fluid simulation is done by 
numerically solving the magnetohydrodynamic (MHD) equations of plasmas 
assuming approximate transport coefficients,while the kinetic simulation con- 
siders a more detailed model involving the particle interactions through elec- 
tromagnetic fields. This is achieved either by solving numerically the plasma 
kinetic equations as the Vlasov or Fokker-Plank equations, or by the newly 
developed Particle Simulation technique, which simply computes the motion 
of a collection of charged particles, interacting with each other and with the 
externally applied fields. 

MHD simulations have generally been applied to large-scale problems directly 
related to the behavior of the experimental devices, whereas, the kinetic simu- 
lations are particularly successful in dealing with the basic problems related to 
fusion physics in which the particle distributions deviate significantly from a 
local Maxwellian distribution, such as in the wave-particle resonance, trapping, 
stochastic heating, etc. The Hybrid Codes developed recently have been the 
bridge between the above two algorithms referred above.In the hybrid codes 
either the fluid or the particle treatment is applied to different components of 
a plasma. 
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1.2 Radio Frequency Wave Heating of Fusion 
Plasma 


Radio frequency wave heating has been proved to be very attractive method 
for bringing the plasma to its ignition temperature in a magnetic fusion device 
from both the physics and the technological stand-points. In a hot magne- 
tised plasma there exists a wide range of energy absorbing mechanisms for 
an electromagnetic wave offered by the plasma dielectric properties. Each 
mechanism is characterized by a perticular range of frequencies of the incident 
electromagetic wave. Depending upon these frequencies, several schemes have 
been proposed to heat the plasma. Each RF heating scheme selectively couples 
energy to a particular species of charged particles (ions or electrons) with a 
desired spatial deposition profile, to certain velocity distribution (Maxwellian 
or non-Maxwellian), and in either the perpendicular or the parallel degree of 
freedom relative to the state magnetic field. 

Besides plasma heating, RF waves offer many other advantages which become 
potentially important considerations in a magnetic fusion device, such as as- 
sisting initial discharge breakdown, driving a plasma current, controlling the 
temperature or current profiles, decreasing unwanted impurities, etc. To date 
there are three major accepted frequency ranges where plasma wave interac- 
tion is found to be strong. The Ion-Cyclotron Range of Frequencies(ICRF^), 
the Lower Hybrid Range of Frequencies(LHRF) and the Electron Cyclotron 
Range of Frequencies (ECRF) . Both the ICRF and LHRF are below lOGHz, 
so generators as well as the transmission systems for such high powerm sys- 
tem are available commercially. The ECRF range for future magnetic fusion 
devices, however, is between 100-200GHz, where high powered sources are still 
to be developed. 


4iere RF refers Range of Frequencies 
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The main advantages of using a RF system in a reactor environment is that 
the wave launchers are thermally shielded from the plasma and the RF sources 
can be placed well behind the radiation shields to avoid direct neutron stream- 
ing. This makes maintaining and operating of high power systems very conve- 
nient. Furthermore, since very efficient RF generators, transmission systems 
and wave couplers exists over most of the RF range, the need for technological 
development in this area is minimized. For RF heating schemes, generators, 
transmission and coupling systems fall in the area of technological studies, 
while the problem of wave absorption is an area of study in fusion physics. 

In the LHRF scheme, high power klystrons are available [if atuan^, (1991)] with 
'^50 kW per tube CW and 200 kW per tube pulsed with typical efficiency of 
about 50-60%. Still higher power, if required (~1 MW) in this frequency range, 
can be supplied by gyrocon which has a very high efficiency (~ 95 %) but much 
lower gain (< 13 dB), so more amplifier stages will be required. Recent gyrocon 
developments have extended the frequency range to 100-200 MHz range which 
is also applicable to higher harmonics of ion cyclotron heating. 

1.2.1 The Lower Hybrid Prequency Range 

In connection with the problem of controlled thermonuclear reactions, RF 
plasma heating has been studied extensively. One frequency range which has 
received a great deal of attention in last two decades is the lower hybrid fre- 
quency range of frequencies or LHRF . The original proposal for the use of this 
method was by Parker. The frequency range has been particularly attractive 
because of following facts: 


• Possibility of plasma heating in large thermonuclear devices. Since fu- 
sion reaction occurs, almost exclusively, among the highest energy ions, 
it is beneficial to have a scheme by which one can produce energetic ions 
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effectively. Ions are involved in wave motion of lower hybrid frequency 
range (as they are massive and lower hybrid resonance is low frequency 
phenomena) and therefore can be heated directly and much more ef- 
fectively than would be the case if higher frequency heating scheme is 
used. 

• Availability of efficient heating power [5^(1965)]. Typical lower hybrid 
frequencies for fusion plasma are of the order of a few GHz, which is 
the highest frequency below which large amounts of power (> IMW) are 
available. On contrary, heating at electron cyclotron frequency, which is 
of the order of 100 GHz, limits one to presently available power of only 
a few kW. Since the heating power required for a tokamak type fusion 
reactors will be of the order of 10® W, this is a very important point to 
be taken into consideration. 

• Lower hybrid waves can be made to heat deep within the plasma as the 
resonance is density and field dependent. 

In last decade much interest has been paid, theoretically, experimentally and 
computationally, to the heating of plasma in d.c. magnetic field at or near the 
lower hybrid resonance [Bekefi{1966)]. It is hoped that sufficient coupling of 
large amounts of RF power to plasma will yield significant increase in plasma 
temperature. In addition, the relatively low frequency of the lower hybrid 
resonance permits large power levels at lower costs. 

Several theoretical models have been proposed to explain the heating process 
at lower hybrid resonance. They are : 

1. Collisional heating, [Bekefi{1966)] 

2. Collisionless absorption [Glagolev{1972a,)] and/or linear mode conversion 
[5^(1965)], 
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3. Nonlinear interactions, especially excitation of parametric modes [Kindel 
et al.(1972)] 

These theories will be discussed in brief to highlight the idea of the problem, 
but before we proceed, a brief discussion of the lower hybrid resonance is 
presented. 

1.2.2 The Lower Hybrid Resonance 

By definition, a lower hybrid resonance is the resonance of the extraordinary 
wave in cold plasrria with the ion-cyclotron frequency. In cold plasma approx- 
imation, the dispersion relation for waves of frequencies above but near the 
lower hybrid frequency has the form [Bekefi{1966)]: 

(3nl -jnl-\-5 = 0. 


The branch , 


n 


2 

X 


7 + ( 7 ^ - 4/36)* 


will be of interest to us, as the choice yields ss A//? as /5 — )■ 0; where 


P = ex, 


and 


A — RL -)- €||Cx, 


where ex, ^xy and ey are the elements of the cold plasma dielectric tensor 
(discussed later), and are defined as 
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where 


e± = 


ClJ/ 


«ll 




1 


(R-L), 


rj2 


R 

L 

ul 


1 f 

^ a;2 \uj + CkO^k 

1 _ y^lS f 

•“0)2 Vo) - ekO^kJ ’ 

A-KnkZle^ 

TTlk 


The solution of ej. = 0 gives 




0)2 - n|r!f 


where Ho is the net plasma frequency (fig + The lowest value of uj that 

satisfies the above equation for fixed values of Qe, and Ho is the lower hybrid 
resonance frequency ujih- In the above case the propagation angle is nearly 
7r/2. Fig. 1.1 shows a schematic dispersion curve as a function of u {ux 
is the refractive index) for cold plasma near the resonance. As the frequency 
approaches near the resonant frequency , the refractive index of the plasma 
increases sharply and, consequently, plasma can be heated near this region, 
since the phase velocity approaches the velocity of the particles. 

As seen from Fig. 1.1, the resonance frequency decreases with increasing plasma 
density. As shown by Stix [^^(IQfiS)], the wave must be slowed down along 
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Figure 1.1; The dispersion curves^ (n‘l vs. ui) near the plasma resonance, ujj, and u)p2 arc 
the resonance points for cold plasma. 

the magnetic field to ensure that increases with increase in density in a 
uniform magnetic field. Stix also showed that for a sufficiently strong deceler- 
ation, the electromagnetic wave in the plasma layer transforms into a plasma 
wave whose longitudinal electric field (with respect to the wave vector) wave 
vector is considerably stronger than its transverse field. Since the refractive 
index of the wave increases with increasing density in the plasma layer, the 
analysis of the wave propagation can be divided into : 

1. The analysis of the ‘Cold’ plasma region where the refractive index is 
not very large, so that the influence of the thermal motion on the wave 
propagation can be neglected; 

2. The analysis in the ‘Hot’ plasma region, where effects of the transforma- 
tion and absorption are important. 


^Fig. 1.1 derived from [Glagolev{1972a)]. 
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In following section, various theories of lower hybrid heating are discussed in 
briefly in order to heighlight the idea behind this work. 


1.3 Theories of Lower Hybrid Heating 

1.3.1 Collisional Heating 

Proposed by Bekefi [5'tia;(1965)] collisions provide the mechanisms for dissipa- 
tion of the RF energy. The propagation vector obtained in the solution of the 
cold plasma dispersion relation will be complex. The best absorption should 
be at the hybrid layer, since the imaginary part of k is largest there. This 
provides the shortest damping length, which is inversely proportional to the 
collisional frequency. 

1.3.2 Collisionless Absorption 

This mechanism was proposed by Stix [5fe(1965)]. The incident electromag- 
netic wave may penetrate directly or excite other modes at the same frequency. 
The wave is then Landau or cyclotron damped, whenever its phase velocity is 
proportional to the particle thermal velocities, which permits the transfer of 
the wave energy from the wave to the particles. In order to study the colli- 
sionless absorption process, it is necessary to consider the dispersion relation 
for plasma waves in a ‘hot’ plasma in a magnetic field in the electrostatic 
approximation. [Chu and Lee(1971), 5tza:(1965)]. 

1.3.3 Nonlinear Processes 

Nonlinear heating processes can occur if the amplitude of the electric fields 
inside the plasma become large enough. The dominant driving mechanism in 
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this interaction is the relative E x B drift motion of the electrons and the ions. 
This drift mechanism can trigger two kinds of processes which leads to growing 
waves in a plasma. They are : 

parametric instabilities [5eA:e^(1966)] which can occur when the driving 
frequency is very close to the lower hybrid resonance. In this case the 
ion acoustic waves with a difference between the driving frequency and 
the lower hybrid frequency can be excited; 

streaming instabilities [Ott and McBride{1972)] whose frequencies are ei- 
ther close to upper hybrid frequency uiuh-, defined as 

uuff = [fig + rig] 

or close to lower hybrid frequency. 

The parametric instability of interest here appears at the lower hybrid reso- 
nance with ion acoustic waves excited [Hook and Bernabai{1971)]. The thresh- 
old RF drift velocity for parametric excitation, which relates the E x B drift 
velocity, Vd = Erf /B o, to the thermal velocity is given by Kindel [Kindel 
et.al{1972)]. As discussed by Ott and McBride [Ott and McBride{1972)], 
plasma heating may primarily be caused by the ” breaking” of the growing ion- 
acoustic waves until the electron temperature is raised such that the threshold 
condition is no longer satisfied. Ion and electron heating was observed by 
Kindel [Kindel et a/. (1972)] in his computer simulation. It was found that the 
relative heating of ions and electron was dependent upon the magnitude of 
the component of the propagation vector along the d.c. magnetic field. For 
k\\/k > the electrons are heated while the ions remain cold. For 

k\\/k ~ the ions and electrons both interact with field and are 

heated simultaneously. 
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1.4 State of the Art 

Computer simulation has proved to be a useful tool in understanding some 
details of the physics related to the overall plasma heating scheme, especially 
those aspects that are difficult to treat analytically. However, to do this eco- 
nomically with the present day computers one must idealize the geometry and 
reduce the scale lengths to manageable proportions. Since simulation of a con- 
finement device is difficult only a basic study has been s presented so far which 
provided a restricted, but detailed view of the problem. The computer models 
presented until now far were used to calculate the self-consistent particle tra- 
jectories of electrons and ions in a magnetized or unmagnetized slab geometry 
to provide a complete picture of the dynamics of these plasma components 
over a short time scale. The purpose was to identify the various factors that 
may also be present in real fusion reactors, and there by suggest the problem 
areas which need further theoretical study. A brief historical overview of what 
has been accomplished with the particle simulation code in RF heating and 
current drive applications has been given in the published work of V. K. Decyk 
[DecyA;(1986)]. Some parts of this publication which are relevant for present 
work are given below for the purpose of quick reference. 

The very first particle simulation specifically for RF heating was by Kindel, 
Okuda and Dawson [Kindel et a^.(1972)]. The geometry was one dimensional 
and an electrostatic algorithm was used. The purpose was to demonstrate by 
theory and by computer simulation that an electric field perpendicular to the 
d.c. magnetic field and oscillating at frequency just above the lower hybrid 
frequency can excite, at extremely low threshold, decay and purely growing 
parametric instabilities which resulted in short wavelength lower hybrid waves 
and ion waves leading to substantial plasma heating. The ExB drift was 
neglected and ions were unmagnetised. This work was later extended to 2D 
[Chu et a/. (1976)]. Subsequent computer models for lower hybrid heating have 
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all been electrostatic, with 2^ dimensions, bounded in one spatial direction and 
periodic in another. The magnetic field direction was chosen so that ExB 
drifts are in ignore able proportions, which means that most of the parametric 
decay processes cannot occur. 

Later models used antennas to excite the lower hybrid waves. The early simu- 
lation with antenna was by Decyk, Morales and Dawson [Decyk et a/. (1979)]. 
The main motive of the work was to understand how to excite lower hybrid 
waves in bounded homogeneous computer models. Vacuum boundary con- 
ditions were used with a specified external surface charge oscillating on one 
boundary. In the early 80’s the interest was diverted to study the ion heating. 
Using the same model Abe and coworkers applied an external charge den- 
sity, specified as an arbitrary function of the periodic coordinates, oscillating 
in front of one of the conducting walls. Plasma was bounded by conducting 
planes. The plasma density profile was inhomogeneous. The major result of the 
simulation was he observation of the ion-cyclotron damping [Abe et a/.(1980)], 
due to the fact that Uo/Cti was rather small (3 or 5), which resulted in high en- 
ergy tail in the ion distribution. The same model but with inversion symmetry 
in one of the bounded direction was used by [Matsuda et a/. (1980)]. 

Later simulations concentrated on the energy absorption mechanism by ions. 
Karney [Kamey{1978), Kamey{1979)] described that stochastic ion heating to 
be the mechanism of ion-absorption which was justified by the examination of 
test particle trajectories for Wo/flj > 8. In order to be able to launch many 
waves whose phase velocities along the magnetic field are large, Nakajima, 
Abe et. al. [Nakajima et a/. (1982)] modified the model to use difiFerent scale 
lengths in parallel and perpendicular directions. In order to avoid numerical 
instability a finer grid spacing and a higher order interpolation scheme in the 
parallel direction was used specially. Inhomogeneous density and temperature 
profiles were used. They also concluded that stochastic absorption is the major 
mechanism for ino-heating. 
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In the mid 80’s the interest in the RF simulations changed from ion-heating 
to current drives. A number of simulations were performed in this area in late 
80’s and early 90’s. As our interest lies only in the heating of plasma without 
the current drive term we will not be going into the details of the current drive. 

Until now almost all the models were based on the electrostatic/electromagnetic 
particle simulation algorithm in 1, l~, 2 or 2| dimensions. These algorithms 
are discussed in detail in a number of published works and books [Decyk et. 
al.(1979)] Almost all models used the vacuum boundary conditions (i.e. there 
is no charge outside the geometry). In most of the models the ExB drifts were 
neglected to isolate the propagation of the lower hybrid waves and to study 
plasma heating without the complications of parametric effects driven by the 
ExB drifts and scattering by the drift waves [Decyk et aZ.(1980)]. The most 
recent models included the ExB and the diamagnetic drifts by tilting the 
axis of magnetic field in the xz plane [Chu et al.(1976)]. All of these models 
were applied to study the parametric heating of plasma by a large amplitude 
pump with frequencies near the lower hybrid and to investigate the effect of 
the parametric instability on the cross-field transport properties of plasma. 

All these models were based on the Cold Plasma Approximation. The ef- 
fects of finite electron temperature (Tg) and ion temperature (Tj) were taken 
into consideration by adding the thermal correction term to the electrostatic, 
cold plasma dispersion relation, which was derived using the Vlasov the- 
ory(discussed in chap. 2). A considerable simplification was made in the ion 
Vlasov analysis while deriving the dispersion relation by neglecting the mag- 
netic field (i.e., treating ions as essentially unmagnetised, Vi x B = 0). This 
is justified in the light of the fact that 1 in the core plasma region. 

Where pi = vm/Di is the ion Larmor radius. In this case the ion gyro orbits 
are practically straight lines as compared to the wavelength Aj. Thus the scope 
of these models was limited in the sense that the region of propagation in which 
the ‘old’ plasma and hydrodynamic approximations are valid becomes smaller 
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as the ion and electron temperature and plasma density increases. The viola- 
tion of the cold plasma and fluid approximation starts near the ‘cold’ plasma 
resonance. 

By using the hydrodynamic approximation it is shown that if the magnetic 
field exceeds a critical value, which is given approximately by the condition 
(me/mi)Qg > the electromagnetic wave should completely transform into 
plasma wave at some critical density in the plasma layer. After being reflected 
from the transformation point, the plasma wave exhibits a stronger decelera- 
tion as the plasma density decreases. As soon as the phase velocity of the waves 
approaches the ion thermal velocity, Cherenkov absorption of the waves should 
occur. This effect should result in heating the ion component of plasma. In this 
case the hydrodynamic approximation is not valid in calculating either the ion 
or the electron absorption. This arises from the absorption between the phase 
velocity of the wave along the magnetic field and the thermal velocity of the 
electrons. When the phase velocity of the waves is close to thermal velocity of 
the particles, their thermal motion should considerably affect the propagation 
of waves in plasma. In all the models discussed above the phase velocities were 
less than almost three times the electron thermal velocity {vthe = y^^e/^e)- 


1.5 Outline of the Work 

Chapter 2 The linear theory of lower hybrid waves is briefly reviewed. The 
dispersion relation with cold plasma approximation is discussed along 
with the warm plasma correction. Finally, plasma heating via hear mode 
conversion is discussed. 

Chapter 3 Some of the important algorithms used in the simulaation are dis- 
cussed. Various protocals are given along with a brief schematic diagram. 

Chapter 4 The input parameters and the plasma model along with the bound- 
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ary condition are discussed. Finally, results of the simulations are dis- 
cussed. 

Chapter 5 Summary and conclusions of this work are presented. 

Appendix A Poission’s solver algorithms, for one as well as two dimensions. 



Chapter 2 


Linear Theory of Lower Hybrid 
Waves in Maxwellian Plasma 


2.1 Introduction 

This chapter briefly reviews the theory of lower hybrid waves in electrostatic 
Maxwellian plasma. As outlined in Sec. 1.5 the theory of lower hybrid heating 
is classified as; (1) collisional, (2) collisionless, and (3) nonlinear processes. 
According to the nonlinear theory, plasma heating takes place because of the 
turbulence that is created due to enhanced level of fluctuations of the internal 
fields of plasma. These fields are created by the E x B drift motions driven 
by a sufficiently strong oscillating electric field, with a frequency close to one 
of the characteristic resonance frequencies of the plasma. Under the action of 
these fluctuations (or turbulence), the particle velocity distribution function 
changes and energy of particles increases, i.e., an anomalously fast transfer of 
the energy from pump field to particles occurs. This nature of the process 
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has attracted the attention of many researchers, resulting in a considerable 
amount of work in this area. On the other hand the advantages offered by the 
linear theory seem to be more attractive. Some of these are discussed in the 
following text. 

Plasma heating in the nonlinear regime, as discussed earlier in Sec 1.5, requires 
a strong pump field. Moreover, for this case, efficient heating can tke place only 
above a certain threshold. In large steady state thermonuclear devices where 
the energy stored is of the order of 10®— 10^ J and where the containment time is 
often large, it is impractical to propose to heat the plasma by nonlinear effects. 
One important reason is the technical problems that are associated with the 
devices which use high-power pulses to create the requred large amplitude 
fields at extremely high frequencies. More promising methods are therefore 
those in which the collisionless absorption can occur with small amplitude 
high frequency fields. Following points must be satisfied in order to achive this 
type of heating; 

1. Penetration of high frequency fields into high density plasma; 

2. Retardation of these waves to velocities such that collisionless (Cherenkov 
or cyclotron) absorption of the waves is possible; 

3. Maintenance of plasma stability during heating. 

The rest of the chapter is organized in the following manner: 

Sec. 2.2 briefly reviews the theory of wave propagation in Cold Plasma. The 
section is concluded with a discussion of electrostatic and electromagnetic 
cold plasma dispersion relation. 

Sec. 2.3 VFam PZasma effects are discussed. The section also summarizes the 
analytic approximation to electrostatic hot plasma dispersion relation for 
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a Maxwellian plasma in the region of lower hybrid resonance, which in 
particular, is associated with the inequality 

... 1 ... n 

>1 or cjo > “t 

a Li 

Sec. 2.4 Plasma heating in the linear regime via mode conversion is discussed. 

The conditions for the propagation and absorption of electromagnetic 
waves in cold and hot regions are reviewed. 

Sec. 2.5 This section concludes the chpater outlining the discussion by various 

authors as [5tix(1965)], [Glagolev{1972a,), Glagolev{l972h)], [Brambilla{1976)] 
and [Chu et a/.(1976)] regarding the mechanism of energy absorption in 
the linear regime. 

2.2 Wave Propagation in Cold Plasma 

2.2.1 Electrostatic Approximation [Bonoli, (1983)] 

For the situation shown in Fig. 2.1 with ambient electron and ion plasma 
density no{x) and static magnetic field Bq = zBg, an electrostatic oscillation 
is setup with frequency ujq and wave number k = xkx. + zk\\. Without loss of 
generality k has been taken to have only components in x and z directions. 

The frequency ujq is assumed to satisfy fig > Wg > where f2| = qkBo/mk is 
the gyrofrequency of the k^^ species. The perturbed electric field JE is related 
to the perturbed electron density Sug and the perturbed ion density Srti by the 
linearized Poisson’s equation [Sono/i(1983)], 

V • 5E = — {6ni — dug ) . (2.1) 

In the electrostatic approximation, the perturbed electric fields satisfies the 
relations V x 5E = 0 and 5E = — V(^(x, t), where <i>('x.,t) is the electrostatic 
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Figure 2.1: Cartesian geometry for slab analysis of wave propagation 

potential. The perturbed densities must also satisfy the linearised continuity 
equations 


dSrie 

dt 


+ V • (no^Ve) 


86 Hi 

dt 


+ V • {rioSvi) 


0 , 

0 , 


( 2 . 2 ) 

(2.3) 


where 6ve and 6vi are respectively the perturbed electron and ion fluid veloc- 
ities and Veo = Vio = 0 has been assumed {‘cold’ plasma approximation), 
The perturbed densities are assumed to have the waveforms 


5n, = 

(2.4) 

Sm = 

(2.5) 


(2.6) 


(2.7) 

Svi = 

(2.8) 
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where ne,ni,^i,Ve,Ve are constant amplitudes. Now Eqs. 2.1, 2.2 and 2.3 
can be rewritten as 


k'^(j) = 

— {Sui - Srio ) , 
eo 

(2.9) 

SUe = 

— k • dVe, 

(2.10) 


UJo 

SUi = 

7^0 1 r 

— k • dVi, 

0^0 

(2.11) 


where = A;j_ + A:|. For the above equations to be valid, the WKB approxi- 
mation must be satisfied, according to which 


k I > 


1 drio 
Un dx 


The perturbed fluid velocities are solved in terms of (/> using the electron and 
ion momentum equations. 


g((5E + 5veXBo), (2.12) 

q (<5E 5vi X Bo) , (2.13) 

The results for and dv,, correct to order {ujo/Qe) and to order (Qj/cJo), 
respectively, are (in the limit ':$> Qf), 


d5v^ 


rrif. 


rrii 


dt 

d6vi 

dt 


dVey 

5Vez 




kL(l> 

Bo ’ 

ik±(j) 

OJQ fcj.0 
fie Bo ’ 
fit k_i(f) 
oJo Bo ’ 


(2.14) 

(2.15) 

(2.16) 
(2.17) 
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Sviy = 0, (2.18) 

Sy.^ = (2.19) 

Uo Bo 

Using Eqs. 2.14 to 2.19 in Eqs. 2.10 and 2.11, and combining the results with 
Eqn. 2.9 yields the electrostatic lower hybrid dispersion relation, 


where 


q kf. 


kl 


kl 



are the dielectric tensor components 


He 


2 \ 1/2 
noq^ \ 

eorrie) ’ 


is the electron plasma frequency and 




eoTrii ) 


is the ion plasma frequency. 


If Eqn. 2.20 is solved for k\^ we get, 

k\ = -A:Je||/ex, 


( 2 . 20 ) 


(2.21) 

(2.22) 


(2.23) 


so that there is a plasma resonance {k± — >• oo) for = 0 and plasma cutoff 
{kj_ — >■ 0) for €|| = 0. Using Eqn. 2.21, the condition ex = 0 can be solved for 
cuq to define the lower hybrid frequency ujih as 

2 _ 2 nf 

^LH - ^ n,/02) 


(2.24) 
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which can also be written in the form first given by Stix [5'tza:(1962)] as 

1 ^ l__ 1 1 

^Ih ~ nf ■ 

The frequency Q is the same frequency referred by Stix as the ”gm 

gyrofrequecy” , which is nothing but the geometric mean of the ion cyclotron 
frequency and the electron cyclotron frequency. As discussed in most of the 
publications the solutions of Eqn. 2.23 for typical central plasma parameters 
yields the useful limit 


(^l/^ll) ~ (He/Wo)^ ~ 


{rrii/me) > 1 , 


for an electrostatic lower hybrid wave. 

2.2.2 Electromagnetic Case [5ono/?,(1983)] 

The electromagnetic dispersion relation is briefed here to highlight the reason 
for using the electrostatic simulation algorithm, and also for the purpose of 
quick reference in case of future development of the code. 

The electrostatic approximation made in above section is not always valid 
for the lower hybrid wave. For example near the plasma cutoff where kx ~ 
0, the parallel wave number (A:||) of the wave may not be large enough to 
satisfy (/ cc/wq) = (^1 + ^|) (c/wq)^ ^ 1 (so that V x JE = 0 may not be 
assumed). The fey of the wave is fixed by the axial (along Bo) dimension of 
the wave guide launching structure or gnW which can, in general, launch waves 
which satisfy kc/cvo « k\\c/uJo > 1 near the plasma cutoff. 

The situation shown in Fig. 2.1 is again considered. However, the starting 
point for the analysis is Maxwell’s equations. 
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V X (5B 

V X (5E 

The perturbed electric field 5E, perturbed current density ^J, and perturbed 
magnetic field ^B, are assumed to satisfy the waveforms, 


5E 

jggik-x— xcdot 

(2.27) 

(5B 


(2.28) 

<5J 


(2.29) 


= eo-^^SE + SJ, (2.25) 

= -pB. (2.26) 


where E,B and J are constant amplitudes and k = xkx + zk\\ without loss of 
generalization. Substituting in the Maxwell’s equations (Eqs. 2.25, 2.26) and 
combining them we can get the wave equation as, 

-kxkX(5E= (^y)^(5E + za;o(5J. (2.30) 

The perturbed current density can be expressed in terms of the fluid velocities 
as 

5J = Uoe ((5vi — (5ve) . (2.31) 

The <5vi and are solved in terms of <5E using the momentum Eqs. 2.12 
and 2.13. Substituting these results in Eqn. 2.31 we can solve 5J, then using 
Eqn. 2.30 we get the wave equation as [jBono/i(1983)] 

kxkX(5E+ (^^^)V-E = 0, (2.32) 


ex 

0 


-It 


xy 


ex 

0 


K = 


0 ■ 
0 

^11 - 


(2.33) 
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where K is the famous dielectric tensor of Stix [5t«a;(1962)], evaluated in the 
limit Og Wq Qf and the components are defined as follows: 


1 

{R + L), 

^xy — 2 ’ 

(2.34) 

R 

III 


(2.35) 

VoJo + (-k^k / 

T, 

= 

k ^0 

( .... ^0 \ 

(2.36) 

Ju 

\UJq ~ CkClk^ 

P 

fc ^0 


(2.37) 


_ A-KUkZle^ 
rrik 


(2.38) 


Using the dielectric tensor equation Eqn. 2.33 the wave equation 2.32 can be 
written in a more compact form by defining 

D = kxkx 1+ (2.39) 
Then the equations of the wave becomes 

D • 5E = 0. (2.40) 


If 9 is to be the angle between Bq = zBq and n, and if we assume n to be 
increasing in the positive x direction, then Eqn. 2.40 becomes 


ex — cos^ 9 


cos 9 sin 9 


' E,- 


ex-n2 

0 


Ey 

cos 6 sin 9 

0 

ey — n? sin^ 9 




The condition for a nontrivial solution is that the determinant of the square 
matrix be zero, thus 
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e± — cos^ 9 —ie^y cos 9 sin 9 

itxy ex — 0 =0. (2.42) 

v? cos 0 sin 0 0 e|| — sin^ 9 

This condition gives the famous dispersion relation, or the equation for the 
wave normal surface as [5'i?x(1962)]: 


pnl --ynl + S = 0, (2.43) 

where j3, 7 and delta are the coefficients given by: 

/5 = ex sin^ 9 + ey cos^ 9, (2.44) 

7 = i?I/sin^ 0 + e||ex(l + cos^ 5), (2.45) 

(5 = e\iRL, (2.46) 


and R, L and P are as given in Eqs. 2.35, 2.36 and 2.37, respectively. 

2.3 Thermal Corrrection to ‘Cold’ Plasma Dis- 
persion Relation 

The thermal corrections to the cold plasma dispersion relations (Eqn. 2.20) 
are derived using the Vlasov theory ([ii'aT7rey(1986)] and [Krall et.al. (1973)]). 
The electrostatic equations are used as plasma effects should only be important 
in the high temperature regions of the plasma where the electron and ion 
densities are high enough to satisfy (fcc/wo)^ 1- As electrostatic oscillations 
are considered Poisson’s equation is used. Assuming wave forms for 0, Sue 
and Srii similar to those in the previous section, Eqn. 2.9 can be rewritten as 
([5'tia;(1965)] and [5ono/z(1983)]); 


1 + Xe + Xt — 0) 


(2.47) 
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k'^4>eo ’ 


(2.48) 


where Xk is the susceptibility function for the species. The perturbed 
electron and ion densities are calculated from the kinetic theory [Krall et.al 
(1973)], 


5n^ = rio J 6fe (v) dv, 

(2.49) 

Srii = no [ Sfi (v) dv, 

(2.50) 


where Sfg and 5fi are respectively the perturbed electron and ion distribution 
functions. 


2.3.1 Electron Vlasov Analysis [Bonoli,{1983)] 

The perturbed electron distribution is calculated from the electron Vlasov 
equation 


dFe , dFe 


-5- (E + V X B.) . 

me 



Fe is expanded as 


(2.51) 


Fe (X, V, t) = feo (v) + 5fe (x, V, t) , (2.52) 



where feo is the unperturbed distribution, Vthe = [T'elmef'^'^ is the electron 
thermal speed, and 5fe -C feo- The equilibrium and the linearized electron 
Vlasov equations are 
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dfeo dfeo 

dt ax 


rrie 


(v X Bo) 


dfeo 

(9v 


= 0 , 


dfe 

dt 


+ V • 


9x 


rUe GY rrie dv 


(2.54) 


(2.55) 


where Eq = 0 has again been used. As the electrons are strongly magnetized, 
the ambient magnetic field is retained. This fact is justified in the light of the 
fact that electron gyro orbit is very small as compared to the perpendicular 
wavelength, {k±Pe)'^ <C 1. Here pe = Ue/Qe is the electron gyro-radius. The 
inequality [k^pef' 1 is easily justified using k\ « A:j| (fle/tuo)^. 

The linearized equation (Eqn. 2.55) is solved by first transforming v = (ui, Vy, Vz) 
to a cylindrical coordinate system [5ono/i( 1983)] (uj., uy, u^) where 


Vx = Ux COS d, 

(2.56) 

Vy = v± sin 

(2.57) 

= vl + vl- 

(2.58) 


(2.59) 


With these transformations the velocity gradient transforms as 


d 

. d sin9 d 

(2.60) 

dvx 

^aux ux de' 

d 

d r d 

(2.61) 

dVy 


d 

d 

(2.62) 

dvz 

5U||’ 


This transforms the unperturbed distribution function feo{v‘^) to a function of 
(ux,U||) because = v1_+vj. feo{v±,v\\) can also be obtained from equilibrium 
Vlasov equation Eqn. 2.54. Recalling that feo now is independent of x and t, 
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and transforming the velocity gradient, Eqn. 2.54 yields 


Q ^/eo _ Q 

de 


Thus, /eo(^Xj^^ll) is independent of 9. The electron orbits are governed by the 

equations of motion 


II 

(2.63) 

= ^v(t)xBo. 

at rrie 

(2.64) 


Scalar multiplication of Eqn. 2.64 with unit vector b in the direction of the 
magnetic field gives dv\\/di! = 0, where, uy = b • v(t‘). Thus, v\\ is constant of 
particle orbit motion. Let Wk stands for total electron kinetic energy 




(2.65) 


where Wn = Imeuji and W±_ = ~mev'j_. Since Wk must be conserved. 


dWk 

dt 


= me 


U|| 


dt 


+ v± 



= 0 . 


( 2 . 66 ) 


The linearized Vlasov equation (Eqn. 2.55) is now solved by integrating along 
the orbits of the electrons (Eqn 2.63) in the unperturbed fields. Eqn. 2.55 can 
now be written as 


dSfe 

dt' 


[x(t'),v(t'),t'] = —6E[x{t),t] 

mp 


dfeo 

dv 


/eo[v(t')]. 


(2.67) 


where the left hand side is now the exact differential. Being constructed from 
the constants of motion, feo is constant along the orbits. Integrating Eqn. 2.67 
in the limits t' = — oo to l! = t, assuming Sfeo = 0 as — )■ — oo yields 
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<5/e[x(i), v(t),t] = -z— / dt f dT[kx cos 9 {t)^^ + (2.68) 

77lg oo J —oo d'^1. 

where r is defined as t = t' —t. Using Eqs. 2.60 to 2.62 to transform k- (5/5v) 
and performing a change of variable in the integral to t = t' — t v/e have 


dfeo , , dfeo 


(5/e[x(t), v(t),t] = -i—<f)[x{t),t] f 

Trip J-c 


[fcx cos 0{r)^ + /C|| ■ (2.69) 

f -- - -- [sin(Qe'?- — 9{t)) sin0(t)] + j(^||U|i — oJo)t\ . 

Using Eqn. 2.69 in Eqs. 2.49 and 2.49 expression for Xe is obtained as 


exp (i [sin(Qe'?' — ^{t)) sin0(t)] + i(^||U|| 


.Ui /•2’r roo poo rO 

-¥/o ‘^Jo J-J^ 

exp fi-^^[sin{QeT — ^(f)) sin^(i}] + *(^||y|| — , 


(2.70) 


where dvxdvydv^ = d9v±dv_idv\\ is used. Eqn. 2.70 can be written in the 
approximate form, as given by Glagolev [(?/a5o/ev(1972b)] 


where, 


2U^ 


(2.71) 


z 

= X-iY, 

(2.72) 

X 

= 2xe~^^ f e‘'^dt, 

Jo 

(2.73) 

r 


(2.74) 
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X = ujo/kwvthe, lo = h{Ze), Zf. = k\Tg/meZt1 Jo is the Bessel function of 
zeroth order. For wo/^n'Ut/ie ^ 1, Xe is approximately 


Q2 ulk'^° 

(2.75) 

(jJq K 

where 

2.3.2 Ion Vlasov Analysis [5onoZi,(1983)] 

The ion Vlcisov equation is 


dFi ^ dFi 


— {E + vxBo)~ = 0. 

TTLi av 


(2.76) 


The perturbed ion distribution is calculated using the above equation. The 
total ion distribution function is expanded as 


Fi{x,v,t) = /io(v) + 5/i(x,v,t), 

, / 1 \3/2 1 f \ 

with fio as the unperturbed ion distribution, vthi = {KTi/niiy^'^, the ion ther- 
mal speed, 5fi <C fio, and v'^ = vl + Vy+ v^. The linearized Vlasov equation is 
then 


(2.77) 

(2.78) 


dSfi , d5fi , q dfio 

___ 4- V • — 1 Ohi ■ — 

ot ax. rrii av 


= 0 , 


(2.79) 
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where the unperturbed electric field Eq = 0. In writing Eqn. 2.79 a consider- 
able simplification is done by neglecting the magnetic field (i.e., treating ions 
as essentially unmagnetized). This is justified in the light of the fact that the 
lower hybrid resonance region is particularly having very long ion gyro orbits 
as compared to the perpendicular wavelength Ax (i.e., {k^piY :§> 1). 


Assuming that 5 fi can be represented by a waveform 




and 


5E = — V</>(x, i) = ik^(x, t), 


Eqn. 2.79 can be solved for (5/i, 


<5/i = 


rrij Wo — k • V 


(2.80) 


(2.81) 


(2.82) 


Using this results in Eqs. 2.48 and 2.50 yields 



/ 


k-(^) 

Wo — k • V 


dv. 


(2.83) 


The approximate equation given by Glagolev [GZo^o/ev(1972b)] is 


Ai = -2 





where Z is as defined in Eqn. 2.72. 


(2.84) 


In this case Xi does not include the absorption at high harmonics of the ion- 
cyclotron frequency which can occur when Wq ^ f2i, for which the Larmor 
radius is large as compared to the wavelength. 


The imaginary part of the function Z in Eqs. 2.71 and 2.84 represent the 
Cherenkov absorption of the "waves by the electrons and ions, respectively. . 
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2.4 Plasma Heating Via Mode Conversion 

When the phase velocities of the particles are larger than their respective 
thermal velocities, the wave propagation equations for infinite plasma are given 
by [Glagolev{1972&)] 

+ = (2.85) 

SH, + [(j, - 1^) - ij] V, = (2.86) 

Eqs. 2.85 and 2.86 characterize the ‘E’ and ‘H’’ waves. For k\i = 0, Eqs. 2.85 
and 2.86 represent the ordinary ‘E’ and extraordinary ’H’ waves, propagating 
in the direction perpendicular to the magnetic field, i.e., in the x direction. In 
case of finite A:||, coupling occurs between the ‘E’ and ‘E’ waves. Assuming 
k\i = kn\i and k_i = kn_i, and solving the dispersion relation (Eqn. 2.41), we 
get 

^li ,2 =p±\/p‘^-q (2.87) 

where p = 7/2/5 and g = 5/p. The (+) and (-) sign correspond to the ‘E’ and 
‘E’ waves, respectively. Thus the ‘E’ wave, also called the quasi extraordinary 
wave, is a backward traveling wave; therefore the waves propagating from the 
external antenna towards the lower hybrid resonance will have a negative phase 
velocity. As seen from the above equation, the ‘E’ wave has a larger refractive 
index than the ‘E’ wave, and it is the ‘E’ wave that must transform into the 
plasma wave. Figs. 2.2 and 2.3 show the characteristic variation of with the 
X coordinate (representing the direction of increasing density) for the wave of 
‘E’ type in the cold and hot plasma regions. The waves can propagate in two 
different ways as seen from the plots, as 

1. When the electromagnetic wave (‘E’) converts smoothly into a plasma 
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wave (transformation point being x = a, Fig. 2.2); 

2. When there is no smooth transformation of ‘E" wave into a plasma wave. 

The condition for the efficient transformation of electromagnetic ‘E” wave into 
a plasma wave is derived by Glagolev in his work [Glagolev{l972a.)] as 

I 112 n 

The inequality 2.88 is also a sufficient condition for the propagation for both 
types of waves. The separation of the waves is observed to increases sharply 
with increase in ny above the value given by the equation 

n||~n^ = v^+^ (2.89) 

The wave damping can thus disappear in the region between cold and hot re- 
gions if a sufficiently high density is present at this layer. In this case, for small 
densities, the longitudinal deceleration will considerably higher as a result of 
which the ‘E" and ‘iif’ waves will likewise separate at the plasma boundary. 
Since the 'E' waves has a higher deceleration along the x axis inside the plasma 
layer, only that wave will be transformed into the plasma wave. Therefore, it is 
important to generate the ‘E’ wave on the plasma surface if they are to be ex- 
cited inside the plasma. Nevertheless, the excitation of ‘if’ waves on the plasma 
surface can be useful for providing the hydrodynamic stability to the plasma 
layer. Since, ‘if’ waves exist only on the surface skin layer it exerts a high 
frequency pressure on the plasma. According to Glagolev [Glagolev{1972&)] 
this high frequency pressure leads to the hydrodynamic plasma stability. 
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2.5 Summary 


From the discussions in Sec. 2.4, it is clear that the transformation of the 
waves is of great importance for the absorption of waves and therefore the 
heating of plasma. According to Stix [5tia:(1965)] the ‘F’ and ‘if’ waves must 
sufficiently be decelerated along the magnetic field in order to penetrate into 
the plasma with direction of the density gradient perpendicular to the direction 
of magnetic field. .Also the square of the propagation vector reduces as density 
increases, and thus the hybrid resonance layer is not accessible from the lower 
density side of the plasma. In order to avoid this difficulty, a component of 
the propagation vector parallel to the magnetic field is required (i.e., > 0). 


If the magnetic field approaches the value given by 

1 


^ > 

9 




me 


(2.90) 


where Ug = and q*’ is the damping coefficient for the species, 

as given by Glagolev [Glagolev{l972a)], the density corresponding to the wave 
transformation conditions increases without limit. Also, if the magnetic field 
is less than the value specified by Eqn. 2.90, transformation of the 'E' wave 
into plasma wave is impossible. In other words it can be stated that, if the 
magnetic field exceeds the critical value corresponding approximately to the 
condition > loq, the electromagnetic ‘F’ wave will, at some definite 

value of density in the plasma wave, completely transform into the plasma 
wave. After reflection from the point of transformation, the deceleration of 
the ion hybrid waves increases with the plasma density. As soon as the phase 
velocity of the waves approaches the ion velocity, Cherenkov absorption of the 
waves by ions occurs. This effect leads to the absorption of the waves and 
heating of the ion component of the plasma. 


When the phase velocity is close to the ion thermal velocity , their thermal 
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motion should considerably affect the propagation of the waves. The com- 
ponents of the dielectric tensor then depend on the plasma temperature and 
phase velocity and are found by means of the kinetic equation. Glagolev 
[Glagolev{1972h)], has solved Eqn. 2.47 using Eqs. 2.75 and 2.84 and obtained 
the approximate equations 


and 



where + ^ 2 e 3'iid (k = species. In this expression 

u) 2 i and u> 2 e represent respectively the Cherenkov ion and electron absorptions. 
The numerical solutions of this equations (refer, [Glagolev{1972h)]) yields the 
following conclusions: 


1. The critical density increases as the magnetic field decreases; 

2. As ion temperature increases the transformation region shifts to the lower 
density side. This implies that increase in ion temperature improves the 
wave transformation; 

3. The plasma density corresponding to the maximum absorption increases 
slightly with decrease in the magnetic field; 

4. The electron absorption increases with the phase velocity of the wave. 
Therefore, the electron absorption is maximum for small plasma densi- 
ties. 
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5. As electron temperature increases the electron absorption decreases. 

6. Electron absorption decreases as increases; i.e., if the electron Larmor 
radius increases when the wavelength becomes smaller compared to the 
Larmor radius [Z^ ^ 1), the electron absorption decreases considerably. 

7. Electron absorption increases with decrease in the phase velocity along 
the magnetic field. 

This process of mbde conversion is also explained by Stix [5^(1965)], Bram- 
billa [Bramhillail^l^)] and Wong [Wong and Tanp(1976)]. The validity of this 
approach is evident from the form of Eqn. 2.43 in which the magnetic field has 
completely disappeared from the ion term. Indeed for 

fii, which is a characteristic of lower hybrid frequency range, one expects the 
magnetic field to play little role on the dielectric properties of the ion fluid. 
When this condition is satisfied the damping length perpendicular to the static 
magnetic field can be estimated from Eqn. 2.43 by substituating j = 6e — iSi. 
Se includes finite temperature corrections for electrons, in particular electron 
Landau damping along the static magnetic field. On the other hand for purely 
perpendicular propagation the hot plasma dispersion relation admits only real 
roots for and no damping for ions is possible. While, as discussed by 
[Glagolev{1972h)], if the effects of the static magnetic field on the ions are 
negligible, an ion Landau damping is expected by free-streaming of ions. It 
is important to note this conflicting point. This contradiction is explained to 
some extent by the concept of phase mixing, such that, for uq ;:§> Hi phase 
mixing appeares for times -C At <C but after, a whole cyclotron 
period the exact phase correlation is restored, and no damping is left over. 


■O' 



Chapter 3 


Overview of the Algorithm and 
Structure of the Program 


3.1 Introduction 

The idea of obtaining more or less a valid plasma physics using a computer 
to trace the charged particles till their fusion was proposed by many people 
as early as 1950’s . The development of computers and considerable research 
on the the numerical techniques has led to the refinement these algorithms . 
The simulation tools are now capable of producing almost all the electric and 
magnetic interactions as they occur in real plasmas. 

A simulation of real laboratory plasma using a particle approach will lead to 
the obvious difficulty of storing the position and velocity vectors for all the 
particles in the computer’s volatile memory for processing . To avoid this the 
simulation algorithm uses a reduced number of particls each with large charge 
and mass keeping the charge to mass rato neraly same for electrons . This 
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results in increasing fluctuations about the mean (say, of density, ^ l/^/n) 
called as the shot noise. Another difficulty is the non-physical interaction of the 
particles with the spatial and temporal grid used for computing the quantities 
such as the charge densities, plasma potential, electric field and forces, which 
also results in increased noise, called the, grid noise. If the noise level rises 
above a limit it can mask out the desired simulation. 

With the present day fast processors and volatile memory (or RAM) available 
in the range of few Mbytes to something around Gbytes one can attempt to 
reduce these noise levels by increasing the number of particles and using a 
finer grid. This can be easily applied in one dimension, can be pushed to 
two and two and ane-half Dimensions, but presently seems impossible in three 
dimensional simulations. 


3.2 Validity of The Algorithm 

An obvious ■ question that may be asked is how a valid physics can be ob- 
tained by simulating a real plasma with a density of about 10^^ using a 

model which uses a few hundred particles? Before answering this question the 
definition of plasma has to be considered . 

Definition of Plasma 

A plasma is a conglomeration of positively and negatively charged particles. 
It is on an average neutral, as the total negative and positive densities are 
same. The plasma is said to be fully ionized if it does not contain any neutral 
species. A fully ionized plasma has two kinds of charge carriers : free electrons 
and positive ions. Sometimes it also contains negative ions and more than one 
kind of positive ions. Ions may be singly charged or multiply charged; they 
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may also be atomic or molecular. 

Any ionized gas can’t be called a plasma, as there is always some small degree 
of ionization in any gas. So a complete definition of plasma can be constructed 
as follows : 

A plasma is a quasineutral gas of charged and neutral particles which exhibits 
collective behavior. 

By quasineutrality we refer to the fundamental characteristic of the behavior of 
plasma to shield out electric potentials that are applied to it. For the purpose 
of definition it can be stated that if L a characteristic dimension the system, 
is much longer than Debey length A© , then whenever a local concentration 
of charge arises or external potentials are introduced into the system these 
are shielded out at a distance short compared with L leaving the bulk of the 
system free of electric potentials. 

The term collective behavior comes into picture as a small density perturbation 
or local concentration of charge produces electric and magnetic fields within 
the system and affects the motion of other particles. 

Physically, electrons and ions are point particles for most of plasma physics. 
Some interparticle interactions occur at short range in short times, produc- 
ing effects at short wavelengths and high frequencies. Fortunately the details 
of such interactions are of little interest for many purposes. The more im- 
portant effects are the collective, long range interactions which produce wave- 
lengths much larger than interparticle spacing and frequencies with periods 
much larger than the time taken by the particles in crossing the grids. As 
we are more interested in the collective behavior of the plasma,at wavelengths 
larger than the Debey length A > Xo- And secondly the models used are 
intended to produce the essence of the plasma without all of its details. Since 
we simulate over limited time and space , small errors are tolerated and use of 
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ion-to-electron mass ratios {rriilme) like 10 or 100 is used. 


3.3 Overview of the Program 



Figure 3.1: A typical time cycle in the program. The time loop is repeated for each 
desecrate time step till a specified time is reached. 


The program discretizes time and proceeds discontinuously with finite time 
steps, thus using digital rather than an analog computation. In each time 
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cycle the fields are calculated about spatial the grid points solving Maxwell’s 
electromagnetic equations, knowing the positions and velocities of the particles. 
These calculated electric and magnetic fields are substituted in the Newton- 
Lorentz equations of motion to calculate the forces on the particles. The forces 
are then used to calculate the velocities in next time step and hence the po- 
sitions. This procedure is repeated at every time step (see Fig: 3.1). The 
time cycle starts at t = 0, with specified initial conditions for positions and 
velocities . The program runs for a specified number of time steps. 


3.4 Particle Weighting 

At any time step the program calculates the fields at the grid points by col- 
lecting the charges around the grid. This procedure of collecting the charges 
from the particle position coordinates is called weighting. The weighting may 
be of first-order, second-order or higher-order depending upon the order of 
the equation used in the interpolation. A number of weighting schemes are 
available for any order. PIC or Particle- In- Cell being the most widely used 
first-order scheme is followed in this program. 

3.4.1 The Particle-In- Cell Method 

Originally developed by Frank Harlow [Harlow{1964:)] and collaborators at Los 
Alamos in 1955 for doing multidimensional compressible fluid computations, 
this scheme gained a wide popularity in the particle simulation area. The basic 
ides of PIC is to combine the best of the Eulerian and Lagrangian features. 
The numerical instability and mass diffusion of the Eulerian methods and the 
cell distortion difficulties of the Lagrangian methods are overcome. The region 
of interest is divided into Eulerian cells for purpose of computing the field 
variables and the material is transported from cell to cell, as in the Lagrangian 
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Area assigned to ’A’ 



Area assigned to ’B’ 


I - I Area assigned to ’ C’ 
I I Area assigned to ’D’ 



(a) 


(b) 


Figure 3.2; PIC weighting used in charge assignment showing the (a) area assigned to the 
various grid points in 2D, and (h) effective particle shape S{x) in ID. 


approach, in the form of disecrate simulation particles. In the compressible 
fluid computations, where this scheme was applied first, the velocity is a single 
valued function of position and therefore is a grid quantity. On the contrary, 
in collisionless plasma computations their is a distribution of particles at every 
point in space, and it become important to treat the microstructure in detail 
in order to obtain important physics. Accordingly, velocity becomes particle 
quantity, and a simulation particle represents a number of real particles with 
the same electric charge, mass, position and velocity. The grid quantities 
then are the summed up charges and electric fields calculated using the above 
procedure . 

The PIC scheme is also called the area weighting or bi-linear weighting due to 
its geometric interpretation. As shown in (Fig: 3.2(a)), the weights are given 
as ; 

{Ax-Xk){Ay - Vk) 

~ AxAy 

XkjAy - yk) 

AxAy 


Pj+Uk 


3.4 Particle Weighting 


45 


Pj,k+1 

Pj+l,k+l 


(Ao: - Xk)yk 
AxAy 

XkVk 
Ax Ay 


where, rrfcandyA, are particle positions relative to the grid point pc is the 
uniform charge density filling the cell {J2jqj/A, where j = 1, • • • , total no. of 
particles in the cell and A = cell area). In order to resolve the details deemed 
necessary and avoid non-physical instabilities the grid size must be kept as 
small as possible, whereas, in order to make the program computationally 
possible or inexpensive the grid size must be large. Thus an optimization is 
required to cater to both the needs. 


In one-dimension the PIC is just a linear interpolation between two adjacent 
grid points and hence the shape factor {S{x)) (see Fig: 3.2(b), produced is 
triangular in shape. 

The forces are calculated in a manner similer to the one used in the assignment 
of charges. Considering the fields on the grid points, a reverse interpolation is 
made and forces are assigned to the particles. 

Various other schemes are also available for calculating shape factors. The 
most commonly used schemes in the particle simulation are : 


• NGP or Nearest Grid Point [Zeroth- order) [Birdsall and Lanydon(1986)], 

• CIC or Cloud-in-Cell (First-order) [Birdsall and Fuss(1969)], 

• Multipole Expansion [First-order) [Dawson et al.(1973)], 

• Energy Conserving [First-order) [iews(1970)], etc. 


each having its own advantages and disadvantages. Our program is designed 
to use the NGP, PIC, CIC and Energy Conserving schemes and prompts the 
user for his choice. 
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3.5 Integration of Field Equations 


After the charges are collected at the grid points the electric and magnetic fields 
are calculated using the Maxwell’s equations. Due to the density perturbation 
of electrons and ions the electric field generated is given by Eqn. 2.1: 


V -E 


g47r 

^0 


{rii - He), 


P 


Co 


(3.1) 


where, p is the charge density. Solving Maxwell’s equations for an electrostatic 
approximation discussed in Sec. 2.2.1 we have: 


V X E = 
E = 


dBo 


-V(f> 


0 , 


so that, 


( 3 . 2 ) 


Using Eqs. 3.1 and 3.2 we obtain the Poisson’s Equation, 

- (3.3) 
^0 

- (in ID) (3.4) 

^0 

^ (in 2D) (3.5) 

Co 

Once the charge densities at all the grid points are known, they become the 
right hand side of the Poisson’s equation, 

■S7^(f){x,y) = -p{x,y) (3.6) 

with the following finite deference form ; 


= - 

^ dx^ 

d'^E^ d'^Ey 

H = - 

dx'^ dy"^ 


— 2A^ij) 


(A(?!>ij+i + A^ij-i 

Ay2 




-Pi,, (3.7) 
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The above set of equations is solved for all using appropriate boundary 
conditions (open, bounded, periodic, etc.). The GS-SOR iterative technique 
is used in 2D case; for ID case the equation is simply of the form : 


and is solved by Gauss elimination using Thomas algorithm. 

Once the plasma potential (j) is known at all the grid points the E’s are obtained 


using Maxwell’s equation. 


-V</.(x), 


with the two point. Center Space finite difference form being; 


(E i = (‘/’i+lj 

^ 2Ax 

^ 2Ay 


(3.10) 

(3.11) 


This completes the integration of the field equations. The model being elec- 
trostatic (5Bo/5t ~ 0), only static magnetic (Bq) fields are to be used in the 
Newton-Lorentz equations of motion. 

3.6 Integration of Equations of Motion 


The charge density at the laboratory coordinate x of a cloud whose center 
is at X is changed from q6{x — x) for a point particle, to qS{x — x) for a 
cloud, where q is the total charge given by qJdxS{x! — x). If Jp and Pp are 
the current and charge densities of a system of point charges located at the 
(x ), then the densities Jc and pc for a system of clouds, whose centers coincide 
with the point particles are [Birdsall and Langdon{l986)]: 


Pc(x, t) 
Jc(x, t) 


dxS{x - x) 


Pp(x',i) 
Jp (x , t) 


(3.12) 
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These clouds, when used in Maxwell’s equation, give the electric (E) and 
magnetic (Bq) fields. The Newton-Lorentz force on one cloud of total change 
q, with (center) position x and velocity v is then [Birdsall and Langdon{1986)] 

F(x,v,t) = qj dx'5(x' -x)- [E(x',t) + vxBo(x',t)] . (3.13) 

In simple form this equation can be written as (which essentially the form 
given by Eqn. 2.13 


dv 



(3.14) 

dx 



(3.15) 


where, 

^ B electric "b F magnetic 

= q{E + vx Bo) 

In our case (see Fig: 3.3), the position is resolved in x and y directions 
and velocity is resolved in v^, Vy and components. In order to include the 
E X Bq drift in the model, so as to study the dynamics of ion component of 
the plasma, the magnetic field is tilted with respect to the 2 axis in the xz- 
plane. In this case the static magnetic field Bo is making an angle d with the 
2 direction. The self-consistent electric field E and wave vector k are along 
x-axis. Then U|| is the parallel velocity of the particles (along Bo) and vx_ is 
the perpendicular velocity (along x direction). At this point it is convenient 
to solve the motion in the parallel and perpendicular directions with respect 
to Bo which leads to the invention of a new coordinate axis x perpendicular 
to Bo and making an angle 0 with x-axis. Due to this third axis this model is 
sometimes called a two and one-half dimensional model. 

In this new coordinate system the Eqs. 3.14 and 3.15 are solved using com- 
monly used integration called leap-frog method. In the leap-frog scheme the 
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Figure 3.3: The two and one-half dimensional simulation model for the lower hybrid 
heating in linear regime. External Eq is along y direction, Bq in the xz plane. Self consistent 
held in the xy plane. 

finite-difference form of the equations is [Birdsall and Langdon{1986)]: 

Vt+At-Vt _ 1^ 



At m 


(3.16) 

(3.17) 


where Vt+At and Xt+At are new values of velocities and position respectively in 
time, as against vt and Xt which are old and therefore known values. 

In the leap-frog integration there is a phase difference of half a time step in 
position and velocity. Thus the velocity is lagging by half a time step with 
respect to the position. Thus in order to adjust the kinetic energy (due to v) 
and the potential energy (due to x) at same time an additional work is to be 
done in shifting the velocities by half a time-step forward using the force F 
calculated at t = 0. 
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With only a static magnetic field Bq (as our model is electrostatic) the term 
q{y X Bo) in the Newton-Lorentz equation is just a rotation of v, i.e., v does 
not change in magnitude. However, qE = xqE^ + yqEy does affect the and 
Vy components of v. Hence a physically reasonable scheme, proposed by Boris 
[j5om(1970)] is used. It is centered in time as follows : 

The equation of motion for the particles to be integrated is 

^ = 1(E + vxBo). 
at m 

The finite-difference form of Eqn. 3.18 using Center-Time scheme is 

Vt+At/2 ~ 'Vt-At/2 q + Vf_At/2 _ '\ 

At = + At 

The vector equation for Vt+At /2 can be solved for three simultaneous scalar 
equations, one for each component. Other than this the the method suggested 
by Boris, which separates the electric and magnetic forces completely is used. 
Accordingly, 

Vt_At/2 = (3.20) 

' m 2 

Vt+Af/2 — H — . (3.21) 

m l 

Substituting these equations in Eqn. 3.18 

~ i- (v+ 4- V-) X Bo (3.22) 

At ^ 

This leaves only the rotation of v. The equations of motion in the {x , y, Bq) 
coordinates are then integrated in three steps: 

1. add half acceleration to Vt_At to calculate v", 

2. rotation of velocity as per Eqn. 3.22 to calculate V*", 

3. again adding an half acceleration to calculate Vt+At- 


(3.18) 

(3.19) 
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This is commonly called as accel-rot-accel scheme. 

Consider the case where Bo is parallel to z axis. In the x — y plane the rotation 
is through an angle 9 where, 

9 qBo At 

tan - = . 

2 m2 

In the electrostatic case Bq is fixed, thus, 

9 /qBoAt\ ( At\ 

2 \ m 2 J \ 2 J 


Let 


tan 


2t 


then using the half angle formulas, we have 

s = — sin 0 = 

c = cos 9 = 

l+t2 

Then as per the Boris scheme the rotation becomes 


l+t2’ 

1-t^ 


c 

s 




' 1 

—s 

c 




1 

1 


In our case when the directions of Bo and v are arbitrary, the method proposed 
by Boris is used. First v~ is incremented to produce v which is perpendicular 
to (v+ + v~) and Bo. Then 

v' = v“ + v“ X t. 

Since the angle between v“ and v' is 9/2, we have 

. 9 qBo At 

t = -btan- = 

2 m 2 

Finally, (v*" + v~) is perpendicular to v' x Bo, so 

v''' = v“ + v' X s, 




A 123284 
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where s is parallel to Bq and its magnitude is determined by the requirement 
I V- p=| 1^: 

2t 


The scalar equations of this vector form in {x,y,Bo) can be written as; 

= Vx{t- At/ 2) + — cos B 
'‘m2 

Vy, = V,(t~/^t/2)+l^Ey + ^^E,„ 

Til Zi TTl JL 


Vj {t + At/2) 

Vy{t + At/2) 

U5„(t + At/2) 


q At ^ 

cv ' + svy^ H —Ex cos 6 

q At 

-SV^> + CVy^ q —Ey + 


At 


m 2 


E, 


ext 


v^^{t-At/2) + ^^ExSm6. 

TTl ju 


As discussed above the complications arising due to the initial conditions for 
X and V being specified at the same time, i.e. at t = 0, and the main loop 
running with v leading x by half the time step, i.e. by At/2, are resolved by 
first rotating v(0) by an angle AB = +VlAt/2 and applying half acceleration 
using —At/2 based on E(0) obtained from x(0). This point is discussed in 
detail in [Birdsall and Langdon{1986)]. 


3.7 Initial Particle Loading 

Starting with prescribed density profiles for position no( x,y) and velocity 
/o(v), i.e. the velocity distribuation function, and assigning initial particle 
position (x,v)i, is termed as initial particle loading in particle simulation lan- 
guage. The formal process is called inversion of the cumulative density. 
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3-7.1 Position Loading 

Assigning particle position . (a:,y)i from the known density is called Position 
loading. As a weakly homogeneous plasma is considered in this model, the 
density distribution is a linear function of x while it is constant in rest of the 
coordinate directions. The linear function must be such that a variation of 
2:1 is achieved in the positive x direction. This function is generated by the 
program using equation 

no(x) = a(1.0 + ar/Li), 

where a is calculated such that if n^o, %o are total number of particles in 
respective coordinate directions, 

Lx 

= n^o- 

x=0 


Thus, a general form of equation can be written as 

Tt ! '^yOi 

Xi = Li/no(x). 

3.7.2 Velocity Loading 

The model considers only the Maxwellian or Gaussian particles where the ve- 
locity distribution in v is Maxwellian (or Gaussian) of the form exp(—v^/2v'ff^). 
Then using the method of inverting Gaussian the velocities are assigned to the 
particles. A thermodynamic equilibrium distribution, which is Maxwellian, 
is of the form given in Fig: 3.4. Around 99% of the particles are in region 
= 3 and thus we place all the particles in the region v = vth- 

Let Rg be a set of uniformly distributed random numbers varying from 0 to 1. 
Then the cumulative distribution function for u =| v | is 
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where, Rn 


— Ri ~ 6.0), and i?j is a random number between 0 and 1. 


3.8 Organization of the Program 

The program is divided into different modules as per their functions; each 

source module is therefore dedicated for one particular operation. The details 

of various modules and a brief schematic flow-diagram is as follows. 

Input The required data is read either from the input file or from the screen 
as specified. The input file is created using a special program module 
input which is run separately. 

Init After the data is read the module init is invoked which does the initial 
position and velocity loading. 

ranf Assigns a Maxellian’s distribution to the particle velocity profile. 

nomi-vx Here the velocities are shifted half a time step back on time scale. 
This module invokes module advance for the purpose. Advance is thus 
shared by norru-vx and accel. Thus the input values assigned to advance 
are important, and accordingly advance shifts velocity either forward in 
time or backward. 

bnd-condn Boundary conditions are applied using this module. This module 
is designed to apply either singly periodic, doubly periodic or open circuit 
boundary conditions. 


set^rho Set_rho is used to collect charges as per the weighting specified. 
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fields The grid quantities are calculated using fields. Fields invokes module 
Poisson-solver to solve the field equations. 

Poisson.solver The Poisson’s equations are solved using this module. For 
2D case it uses a GS-SOR iterative method and for ID Thomas algorithm 
based on the Gauss-elimination is used. 

advance As per the weighting specified this module calculates the force on 
each particles. Then substituting the electric and magnetic fields calcu- 
lated earlier in other 'modules it solves the Newton-Lorentz force equa- 
tions to calculate the velocities in next time step. 

accel Calculates velocities in next time step by invoking advance. 

move Advances the position using the new velocity calculated in module 
accel. 

The program starts by reading the input either from a keyboard or from the 
input file created by module input (see Fig: 3.5). At f = 0 modules init, 
norrn.vx and bnd-condn are executed for each species. Modules seLrho and 
fields are then invoked to calculate initial fields. 

After completing the initial processing at t = 0 the program steps into the time 
loop. At the start of the time loop module accel is executed followed by move. 
Module plots does the job of making a copy of specified variables at specified 
time steps on a magnetic tape. Finally the charges are recollected knowing the 
new positions and velocities. The time iteration ends with calculating fields 
at grid points; the program then starts the next time iteration beginning with 
accel. This process is repeated for a specified number of time cycles. After 
coming out of the time loop, the program ends making history plots for kinetic 


energy. 
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Figure 3.5; A schematic diagram showing a typical run. The names in the boxes represents 
the actual hlenames of various modules used in the program. 
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3.9 About the Program 

The program is developed in C++ to take the advantage of the Object Ori- 
ented Programming or OOP’s facilities provided by it. The compilation was 
done using GNU C++ compiler g-h-h. Error handling was done with LINT 
package on Digital Unix version V3.2C (rev. US). For the purpose of fast 
compiling the use of archives was made. The program modules are listed in 
the archive file named mainlib.a. The code processing was done via SCCS 
(Source Code Control System) to keep the record of major changes in the 
source files. Debugging was done in GNU debugger gdb. Finally the platform 
used was Digital Unix version V3.2C on DEC-Alpha dual processor machine. 

At various {joints in the code the concept of dynamic memory allocation is 
used so as to save the floating memory (RAM) consumption. This increases 
the speed considcu’ably and also improves the memory mapping performance 
of the {jrocessor thereby increasing the floating point accuracy. The output 
files generated for various variables are according to the specification of various 
gra{)hic. |)a(+ages, so that tliey can be directly loaded on any of them. The 
program offers the facility of producing directly the input file for Gnuplot, a 
generally used grai^hic package on UNIX operating systems, for plotting the 
phase-space plots in one or two dimensions in space and any one dimension in 
phase. 


•o 



Chapter 4 


Discussion of Results 


4.1 Computer Model 

An two and one-half dimensional (v^, Vy, v^, x, y) electrostatic particle sim- 
ulation code has been used. The field quantities and the particle positions are 
two dimensional, while all the three velocity components are retained in the 
particle orbits. Doubly periodic boundary conditions were used. The fields are 
periodic in the direction of the external applied (static) magnetic field(i) and 
bounded in the perpendicular (x) direction. For the purpose of this simulation, 
finite-sized (Maxwellian) particles are used. 

For the simulation results presented here, vacuum boundary conditions are 
imposed in the x direction at rr = 0, i.e., no charge exists outside the bound- 
ary. The model is made periodic in rest of the boundaries. The physical 
consequence of these boundary conditions is that, the charges which cross one 
boundary (say, at y = Ly) are represented at the previous boundary (i.e., at 
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y = 0) as image charges with same parameters as velocity, etc. Thus, the 
plasma model in this ciise represents a slice of an infinite plasma. Also with 
such boundary conditions we can avoid the effects of surface waves and sheaths, 
and therefore can isolatci tlie lower hybrid waves. 

The external electric fields are applied as oscillating charges at one of the 
boundary. In this model we have specified these sources at a; = 0 and due to 
the periodicity they appear as image sources at a; = L^,. The static external 
magnetic: field is si)ecified at an arbitrary angle 0 with respect to the i as shown 
in Fig. 3.3. The* full dynamics of the particles are used and there is no guiding 
centre* approximation made in this code. Finite-sized particles are followed 
with their self consistent fields plus the externally applied driving fields. 

The initial ifiasma. density is taken uniform in y and weakly inhomogeneous 
(with a (h'usity variation of 2:1) in x direction having the maximum density 
at X = Lx- Dillerent, initial tem])arature profiles are used for both the species 
and are uniform in space*. These are specified as temperature ratios. The 
initial ('l(*ct,ron tli(*rmal v(*locity was choosen to be Vthe = 0.75(ine, where 8 is 
the finite grid spacing used in the simulation, and He is the average plasma 
freciuency. Tlie dimensionless, average Debye length is then 0.75, where 
A/j^ is defined as A«,, = •(;{/,« /^Tlg. The ratio of electron cyclotron frequency 
to average plasma freciuency (Dg/Ile) of 1.5 is used. Except when otherwise 
specified, ion to electron mass ratio {mifrue) of 100.0 is used throughout. A 
time step of 0.2511,7* is used to filter out the high frequency phenomena and 
isolat,e the low freciuency, lower hybrid, physics. The reflection of the particles 
implies that thcjre is no mechanism for net energy loss in the system. Number 
of grid points used is Lx/ 5 = Ly/5 = 64. 

The external source in this simulation is determined by specifying an external 
surface charge density oscillating with a specified frequency uq on one of the 
boundary. Because the magnitude of the charge density is fixed externally, such 
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a source corrcspoiKls to a constant current driver. The phase of oscillation is 
such tliat tlie ii('l(ls gcuierated by these charges are zero initially. This is also 
called as a quasi- adiabaUc turn on. 


4.2 Organization of the Results 

The results are organized in the form of different sets, each set containing the 

results collected for different parameter values as given below : 

Set I K(!<(i)iug the temperature profile same results are presented varying the 
ion to ehictron mass ratio to see the effect of using different mass ratio 
on the natnr(! of the curve. 

Set II R(‘sults an' pn'sented lor different exciter frequencies keeping the ion 
to eh'ctron mass ratio and electron to ion temperature ratio constant. 

Set III 'Id se(’ tlu' eflV'cU. of tcunperature on heating of various species, the 
mass ratio is k(q)t constant while the temperature ratio is varied. 

Set IV Rx^sults ar(' coll(Kd,(!d for various values of magnetic fields. 

Set V Results arc; presented regarding the spatial variation of kinetic energy 
for both the plasma components. 


4.3 Simulation Results 

As tlie aim of this work is to study the mechanism of the heating phenomena 
in plasmas, the results are discussed on the basis of the time evolution of ki- 
netic energi(is of the species. Instead of considering the total kinetic energy, 
the kinetic (uuu'gy calculated using the most significant velocity component 
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(i.e., -uii or Vs.) of the species is considered. Thus, for the ion case the perpen- 
dicular kinetic energy {Wj, = \{mvl)i) is considered, while for the electron 
case the parallel kinetic energy (W|| = |(mi;|)e) is considered. Energies in 
the simulation are normalized to initial total electron thermal energy. Thus 
[Decyk et aL(1980)] 





1 

ZNpKTe^" 


E(i)j 


1 




If r is the nondimensional time then r = tUg. Thus, the abscissae physically 
represent the number of plasma cycles at electron plasma frequency, while the 
ordinates represent the kinetic energy calculated using the msot significant 
veocity component for the respective species, as a percentage of the total initial 
kinetic energy. 


4.3.1 Set I 

As we simulate over a limited time and space small errors are tolerated and 
ion-to-electron mass ratios (mj/me) like 10.0 or 100.0 are used. Even though, 
to see the variation of different mass ratios on the nature of heating curve, 
following results were collected. Keeping the temperature ratio same, Tg/Tj = 
16.0, mass ratios of rrii/me = 10.0 and rrii/me = 100.0 were used. Other 
parameters were, k\\/k — 2.5(me/mi)’^/^, wo/Ile = 0.15, ivl/uig = 1.0. fig = 
>^Dr = 2-0, El/AirriTg = 1/200.0. Number of particles per species used 
in the simulation were 2 x 128 x 128 on a 64-cell grid. 

Fig. 4.1 to B^ig. 4.4 show the results of simulation using a weak pump where the 
pump energy is small as compared with the initial electron thermal energy (a = 
150, where alpha is the ratio of initial electron thermal energy to pump field 
amplitude, both measured in energy units). For the electron case (Fig. 4.1 and 
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Fig. 4.3), no significant variation in energy is observed for the different mass 
ratios used. For ions (Fig. 4.2 and Fig. 4.4), even though, is observed to be 
lower tlian tlie i)revious case, the nature of variation of the ion perpendicular 
kinetic energy is similar for the two mass ratios. The difference is because of 
the different mass ratios used. 

4.3.2 Set II 


This set contains results obtained using different exciter frequencies. Keeping 
the mass and temperature ratios constant results are obtained for different 
exciter frecjucncicis. The parameter set is as in set I with Elli-KuTe = 150.0, 
= 2.0, 'nii/rn,, = 100.0 and T^/Tj = 16.0. 

It is observed that^ when uq > Ug, the ion absorption is hampered (Fig. 4.6 
and Fig. 4.8), while the electron absorption increases (Fig. 4.5 and Fig. 4.7). 

4.3.3 Set III 

The initial temperatur(i profile does have an impact on the absorption of the 
electromagnetic energy by a species. For higher temperature ratios it is ob- 
served that the electron absorption is weakened (Fig. 4.11). On the other hand, 
no significant difference is observed for the absorption by the ion component 
of the plasma (Fig. 4.12). 

The parameter .set used was, k\\/k = 2.5{me/miy^‘^, too/^e = 0-15, He = 
Oe, = 1.0, A/j, = 2.0, El/iirnTe = 1/150. Number of particles per 

species used in the simulation were 2 x 128 x 128 on a 64-cell grid. 
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4.3.4 Set IV 

As discussed in Sec. 2.4, an electromagnetic wave can propagate in two dif- 
ferent ways inside the plasma: (1) when there is a smooth transformation of 
the electromagnetic waves into the plsama waves, and (2) when the waves 
do not transform. This mode conversion process is of great importance in 
the absorption, and therefore, heating of the plasma. Glagolev, in his work 
[G/ago/et;(1972a)] has derived the condition for the efficient transformation of 
the electromagnetic wave into the plasma wave depending upon the external 
static magnetic field. This is given by Eqn. 2.90. Thus, in this set of results 
we hav(' tried to verily this inequality. The aim was not to find the critical 
magnetic fit'ld, neither to verify the validity of Eqn. 2.90, but to see the nature 
of the heat.ing curv(' in case when the magnetic field is very small (or, in fact, 
negligil>le for t.h<^ simulation purpose). 

Fig. 4.14 and Fig. 4.13 shows the variation of total electron and ion kinetic 
energies with time. Major observation here is the absence of the resonance 
heating of eitlier siiecies. Tlie total heating observed can be explained to be the 
heating which takers place under the action of pump fields and in accordance 
with the collisional Jonh' heating. Thus, these observations state that the 
transformation of the electromagnetic waves into the plasma has not taken 
place for tins case wlien the. external static magnetic field Bq is very small. 

The parameter set used was, k\\/k = 2.5(Tne/mi)^/^, wo/fle = 0.15, He = 
fie, - 1-0, X,), - 2.0, E'^/AimTe = 1/150. Number of particles per 

species used in the simulation were 2 x 128 x 128 on a 64-cell grid. 

4.3.5 Set V 


In this set of results w(! discuss the spatial dependence of the ion perpen- 
dicular kinetic energy and electron parallel kinetic energy. As discussed in 
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Sec. 2.4, the irradiated electromagnetic radiation propagates inside the plasma 
by first splitting into two waves; the ordinary ‘E’ and quasi-extraordinary ‘i?’ 
wav(?s. Th(i ‘E’ wave liaving a larger refractive index than the ‘if’ waves pene- 
trates the phisma and transform into plasma waves. As discussed by Glagolev 
[Glagolev{1972ix)], this separation occurs near the plasma boundary. The pres- 
ence of ‘f/’ waves on the plasma surface skin-layer leads to considerable heating 
of electrons at this point. 

As seen from Fig. 4.18, the energy gain for electrons is inhomogeneous and 
localized in spacx^ within a few electren Larmor diameters from the source. 
Fig. 4.19 and Fig. 4.20 show the variation in the energy gain with time. It can 
therefore Ix^ (a)nclud(Kl that the electron absorption increases with increase in 
the pluus(' velocity of tlu' waves; and the electron absorption is maximum for 
small plasma dimsitit's. 

Lik(! th(' (4(H4,ron cas(^, the ion spatial kinetic energy distribution is inhomoge- 
neous in sparse. But as discussed by many authors [Decyk et a/. (1979)], surface 
ion heating is not observed; instead, the ion heating was localized almost at 
the centre of the plasma. This distribution is observed to be varying with time. 
After almost six plasma cycles the distribution profile is uniform in space. 

At the end the phase-space plots for parallel and perpendicular velocity com- 
ponents of electrons and ions respectively are presented (Fig. 4.21 to Fig. 4.29). 
The the cones obstu'ved for both the species are most propbably the resonance 
cones as predect(!d by tin; cold plasma theory. Fig. 4.21 and Fig. 4.23 show 
the sheatli region arround the source. 
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Figure 4.13: (set IV) Time evolution of the electron parallel kinetic energy (W|p for a 
small external static magnetic field Bo- 



Figure 4.14: (set IV) Time evolution of the ion perpendicular kinetic energy (W]^) for a 
small external static magnetic held Bq. 
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4.3 Simulation Results 



Figure 4.20: (set V)Spatial dependence of electron kinetic energy parallel to magnitic Held 
at 400^^ time iteration (source at x = 0). 
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Figure 4.21: Phase-space plots (isometric projection) for parallel velocity component of 
electrons at tHe = 400 



Figure 4.22: Phase-space plots (orthographic projecion) for parallel velocity component 
of electrons at tHe = 400 
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Chapter 5 


Conclusions 


With refrence to the results discussed in the previous chapter, following con- 
clusions can be drawn: 

1. Transformation of forward travelling wave into a plasma wave is nec- 
essary for heating the plasma. Continuous transformation of the electro- 
magnetic ‘ E' wave into a plasma wave occurs only when the magnetic 
field exceeds a critical value. 

2. Surface electron heating is observed which implies that the separation of 
the electromagnetic wave takes place at the lower density, as predicted. 
The ‘ E' wave having higher refractive index penetrates into the plasma, 
while, the ‘ff’ wave is absorbed at the surface. Electrons having higher 
thermal speeds are responsible for damping of this wave, which results 
in electron surface heating. 

3. As the resonant on heating is observed to take place after some plasma 
periods (i.e., almost after tUe > 200), it can be concluded that the 
heating is because of the resonance that takes place between the reflected 
electromagnetic wave and ions. This implies that a sufficient deceleration 



of the electromagnetic wave is required so that the wave phase velocity 
almost matches the ion thermal speed. 

4. It is ol)S(U'ved that the transformation of lower hybrid waves to plasma 
waves is temperature dependent. The efficiency of the transformation is 
observed to decrease with increase in electron temperature. 

This method of heating plasma near the plasma surface guarantees an efficient 
absorption of the plasma wave if an electromagnetic wave (of the electrostatic 
type) is generated at the plasma surface and sufficiently decelerated along the 
magnetic field, provided the magnetic field is more than the limit proposed 
by Glagolev [Glagolev{\^72^)], and the exciter frequency is almost equal to 
the ion plasma frequency. Since the transformation condition for the waves 
depend only weakly on the plasma temperature (i.e., the ion temperature), 
this method dose not place any practical limitation on the ion temperature. 

The electron heating takes place at the surface while the ion heating is ob- 
served at the inner regions of the plasma. The presence of the outer hot 
electron sheath shields out the hot ions from impurities coming from the walls. 
Moreover, the existance of ‘iT waves on the surface skin-layer of plasma exerts 
a high frequency pressure on the plasma. This can, in particular, be useful in 
providing the required hydrodynamic stability of the plasma layer. 



Appendix A 

Algorithms Used in the Code to 
Solve the Poisson’s Equations 


Poisson’s equation in the following form is used in the code: 


vV = 


€o 


(A.l) 


One-dimensional case Eqn. A.l can be written as 


dx^ 


P_ 

5 

Co 


(A.2) 


then, for any grid point ‘z’ the discretization (center space) is given by 

4>i+\ + 4>i-\ ~ 




A 

eo’ 


(A.3) 


Two-dimensional case Eqn. A.l can be written as 




£_ 

^0 


(A.4) 


then, for any grid point (i, j) the discretization (center space) is given by 
4’i+l,j ~b ~b (pij—i ‘^4*1, j Pij 


(A.5) 
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For the two cases considered above two different schemes are used. The one- 
dimensional case is solved by Gauss elimination using Thomas algorithm, while 
the two-dimensional case is solved using the Gauss-Seidel iterative technique. 


A.l Thomas Algorithm 


The equation 


cpe 

dx^ 


r, 


(A.6) 


which is of the form given by Eqn. A. 2, is discretized by the center space 
approximation. 


Qt+i + Qi-i ~ 2©i _ 
Ax^ 


this discretization can also be written as 

1 


Aa;2 Aa:2 ^ Ax^ 


(A.7) 


(A.8) 


The above equation cannot be applied for the boundary and is solvable only 
for the interior points. Then, if following is the boundary condition 


di0i = n, 

the Eqn. A.8, written for all the grid points z = 1, . . . ,N — 1,N, gives a set 
of equations which can be represented in the matrix form as 


"di 

0 

0 

0 • 

• •o' 


■ 01 ■ 


’ ri 

b2 

<^2 
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0 • 

■• 0 
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r2 

0 
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da 

03 • 

•• 0 


©3 
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ra 

_ 0 

0 

0 

0 • 

■ • _ 


1 

■■ ® 


. . 


where a, b and d represent above-diagonal, below-diagonal and the diagonal 
elements of the matrix, respectively. This type of matrix is called a tridiag- 
onal matrix. The algorithm used requires 0{N) folating point operations as 
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compared to 0(“A^'*) operations generally required in other algorithms. Thus 
this technique is computationally very efficient. 

The Thomas algorithm is divided into two steps: 

1. Forward elimination. 

2. Back substituation. 

Forward elimination The following set of steps are repeated from « = 1 to 
i = N -1 

factor = hi+iM (A-IO) 

di+i = di+i - factor X Oi (A. 11) 

Ti+i = n+i - factor X ri (A. 12) 

Back substituation First the value of 0 at i is calculated as 

©V —fNidn 

then the following algebraic operation is repeated in a loop running from 
z = A" - 1 down to i = 1 to calculate the values of rest of the 0’s; 

Qi = {ri - Oi X Qi+i)/di. (A.13) 

Solving by the above technique, the values 0 at all the grid points can be 
obtained, provided the initial condition is defined. 

A. 2 Gauss-Seidel Iterative Technique 

Gauss-Seidel iterative technique, also called the GS-SOR (SOR stands for Su- 
cessive Over Relaxation), which begins with a guess solutions and corrects the 
guess until the system of equations are solved to a sufficient accuracy. 
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Advantage : It can use the squarity of matrix efficiently (does not generate 
unecessay work fis in Gauss elimination). 

Disadvantage : Convergence is not guaranteed for all matrices and the rate of 
convergence depends on the given matrix; operation count is uncertain. 


Consider the following set of equations 


aiiTi+ 012^2 l-ainXn = ^1, 


(A.14) 


+ <^n2^2 • * • + CLnn^n ^n* 


Make an initial guess for all xi: 


X 




? *^71 • 


(A,15) 


Then the new values of Xi are calculated as 


= x[ + rifaii, 


(A.16) 


where 


Z— 1 71 

bi — aijXj — aijXj 

j=l j=i+l 


Vi is called the residual of the equation evaluated on the newest values. 


The SOR technique involves the acceleration parameter ‘w’ such that 

= xf + ujri/aii. 

I < iv < 2 implies over relaxation and, 0 < cj < 1 means under relaxation. 
Over relaxation is used to accelerate the solution and under relaxation is used 
when convergence is not achieved by Gauss-Seidel. The algorithm of this 
technique is as follows: 


1. assume an initial geuss for all XiS, 
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2. for each aq, i = 1 to N do 

Tj = bi — dijXj 
Xi = Xi + Ti/au 

3. repeate the above steps till convergence is achieved. 

Convergence is given by | RMS |< error, where 

RMS == 

Convergence Criterion 

Convergence is gaiiranteed only if the matrix has a diagonal dominance, defined 
as 



1. for every row 

I Ojj 1^ ^ ^ I dij 

j=i 

and j 7 ^ i 

2. for at least one row 

71 

I dfji I > ^ ' I dij 

j=l 

and j ^ i. 


These conditions are sufficient conditions for convergence. 
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